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1 INTRODUCTION 



It is now known that the hydrogen in the intergalactic 
medium (IGM), which became neutral at z ~ 1000 (Pee- 
bles 1968; Zeldovich, Kurt & Sunyaev 1968), must have 
been reionized somewhere between redshifts 6 and 30, the 
lower limit coming from the lack of a Gunn-Peterson trough 
in quasar spectra at that redshift (e.g. Fan et al. 2000), and 
the upper limit from the bound on the optical depth to the 
last scattering surface measured from the cosmic microwave 
background (CMB; Netterfield et al. 2001). In fact, very re- 
cent results (Djorgovski et al. 2001; Becker et al. 2001) sug- 
gest reionization very close to the lower limit of this range. If 
there are large populations of galaxies or quasars at high red- 
shifts, as is predicted by current structure formation models 
(e.g. Benson et al. 2001a) and as confirmed up to redshifts 
6 observationally (Fan et al. 2000; Stern et al. 2000), 
then reionization is most likely to have occurred through 
photoionization, as both galaxies and quasars emit copi- 
ous quantities of ionizing photons (e.g. Couchman & Rees 
1986). Several models of reionization have been developed 
in recent years (Haiman & Loeb 1996; Gnedin & Ostriker 
1997; Chiu & Ostriker 2000; Valageas & Silk 1999; Cia- 



ABSTRACT 

We develop a coupled model for the evolution of the global properties of the intergalac- 
tic medium (IGM) and the formation of galaxies, in the presence of a photoionizing 
background due to stars and quasars. We use this model to predict the thermody- 
namic history of the IGM when photoionized by galaxies forming in a cold dark matter 
(CDM) universe. The evolution of the galaxies is calculated using a semi-analytical 
model, including a detailed treatment of the effects of tidal stripping and dynamical 
friction on satellite galaxies orbiting inside larger dark matter halos. We include in the 
model the negative feedback on galaxy formation from the photoionizing background. 
Photoionization inhibits galaxy formation in low-mass dark matter halos in two ways: 
(i) heating of the IGM and inhibition of the collapse of gas into dark halos by the 
IGM pressure, and (ii) reduction in the radiative cooling of gas within halos. The 
result of our method is a self-consistent model of galaxy formation and the IGM. The 
IGM is reheated twice (during reionization of Hi and Heii) , and we find that the star 
formation rate per unit volume is slightly suppressed after each episode of reheating. 
We find that galaxies brighter than L^, are mostly unaffected by reionization, while 
the abundance of faint galaxies is significantly reduced, leading to present-day galaxy 
luminosity functions with shallow faint end slopes, in good agreement with recent ob- 
servational data. Reionization also affects other properties of these faint galaxies, in a 
readily understandable way. 
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rdi et al. 2000; Gnedin 2000a; Miralda-Escude, Haenhelt & 
Rees 2000; Benson et al. 2001a), many reaching the conclu- 
sion that reionization occurred between z «7-12, although 
large systematic uncertainties remain due to uncertainties in 
the efficiency of galaxy formation, in the fraction of ionizing 
photons that escape a galaxy and in the density distribu- 
tion of ionized gas in the IGM (see, for example, Benson 
et al. 2001a). If this picture of reionization is correct then 
it is clear that the thermodynamic history of the IGM is 
determined by the formation and evolution of galaxies and 
quasars. 



The photoionizing background responsible for reioniz- 
ing the IGM may also act, directly and indirectly, to in- 
hibit galaxy formation, as was first pointed out by Doroshke- 
vich, Zeldovich & Novikov (1967), and first investigated in 
the context of CDM models by Couchman & Rees (1986). 
Galaxies are thought to form by a two-stage collapse pro- 
cess, in which gas first collapses into dark matter halos along 
with the dark matter itself, and then collapses relative to 
the dark matter within halos if it is able to cool radiatively 
to below the halo virial temperature, thus losing its pres- 
sure support. The second stage of the collapse is necessary 
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in order to increase the gas density to the point where it 
becomes self-gravitating relative to the dark matter, which 
is believed to be a necessary condition for the gas to be 
able to fragment to form stars. In the presence of an ion- 
izing background, both stages of this collapse process are 
inhibited, particularly for low mass halos. Firstly, the ion- 
izing background heats the IGM to temperatures of around 
lO^K, and the resulting thermal pressure of the gas then 
prevents it from collapsing into low mass halos along with 
the dark matter. Secondly, the ionizing background reduces 
the rate of radiative cooling of gas inside halos, mainly by 
reducing the abundance of neutral atoms which can be col- 
lisionally excited. Both of these mechanisms will strongly 
inhibit galaxy formation in halos with virial temperatures 
less than ~ lO^K, and so may have important effects on the 
faint end of the galaxy luminosity function and also on the 
properties of the dwarf satellite galaxies of the Milky Way 
and other galaxies. 

There have been many studies of the effects of an ion- 
izing background on galaxy formation, both analytical (e.g. 
Efstathiou 1992; Babul & Rees 1992; Chiba & Nath 1994; 
Thoul & Weinberg 1996; Kepner et al. 1997; Nagashima, 
Gouda & Sugiura 1999) and using numerical simulations 
(e.g. Vedel et al. 1994; Quinn et al. 1996; Weinberg et 
al. 1997; Navarro & Steinmetz 1997), but in most of these 
the ionizing background was simply taken as an external 
input. A few studies have investigated the more difficult 
self-consistent problem, relating the ionizing background at 
any redshift to the fraction of baryons which had previously 
collapsed to form galaxies, and at the same time including 
the effect of the ionizing background in inhibiting further 
galaxy formation (e.g. Shapiro et al. 1994; Gnedin & Os- 
triker 1997; Valageas & Silk 1999). The analytical studies 
have used a wide variety of approaches and approximations, 
but have generally modelled galaxy formation and the ef- 
fects of photoionization only in a very simplified or partial 
way (e.g. for photoionization, either considering only the 
supression of collapse into dark halos, or the suppression of 
cooling within dark halos). On the other hand, the numeri- 
cal studies were limited in the predictions they could make 
about properties of the present-day galaxy population by 
the range of physics included and by their limited dynami- 
cal range. In the present paper, we present a new model for 
the coupled evolution of the IGM, the ionizing background 
and galaxies, based on a semi-analytical model of galaxy 
formation, enabling us to determine in much more detail 
than in previous studies the effects of photoionization on 
observable galaxy properties. Compared to previous analyt- 
ical studies (in particular Valageas & Silk 1999), the main 
improvements are that we have a much more detailed model 
for galcixy formation through hierarchical clustering, includ- 
ing many different processes, and a more accurate model for 
how photoionization suppresses galaxy formation through 
the two mechanisms described above. In particular, the sup- 
pression of gas collapse into dark matter halos due to the 
IGM pressure is modelled based on the latest results from 
gas-dynamical simulations. 

Our starting point is the semi-analytic model of galaxy 
formation described by Cole et al. (2000), which attempts to 
model the galaxy formation process ab initio, in the frame- 
work of structure formation through hierarchical clustering. 
We then modify this to include the new physics we are inter- 



ested in here. We develop a model for the evolution of the 
thermodynamic properties of the IGM in the presence of 
the ionizing radiation background produced by galaxies and 
quasars, the former predicted by the semi-analytic model, 
and the latter based on observational data. We are then able 
to predict the mean temperature of the IGM and the spec- 
trum of the ionizing background as functions of cosmic time. 
We adapt the Cole et al. model to determine the mass of gas 
able to accrete onto each dark matter halo from the heated 
IGM, and to include the effects of heating by the ionizing 
background. Finally, we include a more detailed treatment of 
the dynamical evolution of satellites orbiting within larger 
dark matter halos, including the effects of tidal stripping. 
The approach of investigating the effects of photoionization 
on galaxy formation by using a semi-analytic model was 
previously used by Nagashima, Gouda & Sugiura (1999), 
but they considered only the heating of gas in halos by the 
UV background, and so our current work represents a more 
thorough treatment of the problem, as well as being based 
on a much-improved galaxy formation model. 

There are two parts to this paper. Firstly, we describe 
how the physics of the IGM/galaxy interaction may be mod- 
elled in a simple way. Secondly, we present results from our 
model, focussing on the evolution of the IGM and ionizing 
background and the properties of local galaxy population. 
We briefly comment on how high redshift galaxies are af- 
fected. In a companion paper (Benson et al. 2QQlb), we 
will explore in detail the consequences of our model for the 
population of satellite galaxies seen in the Local Group. 
The remainder of this paper is organized as follows. In 
and §^ we describe in detail our model for the evolu- 
tion of the IGM and galaxy formation. In §^ we present re- 
sults from this model for the evolution of the IGM and the 
population of galaxies at the present day in the currently 
favoured ACDM model of structure formation. Finally, in 
we present our conclusions. 



2 MODEL OF PHOIONIZATION AND IGM 
EVOLUTION 

We use the semi-analytic model of galaxy formation devel- 
oped by Cole et al. (2000) to determine the properties 
of galaxies in the Universe. The model includes formation 
and merging of dark matter halos, shock-heating and radia- 
tive cooling of gas within halos, collapse of cold gas to form 
galaxy disks, star formation from the cold gas, galaxy merg- 
ers within common dark matter halos leading to formation 
of galaxy spheroids, chemical enrichment, and the luminos- 
ity evolution of stellar populations. The fiducial model of 
Cole et al. (2000) (for which D,o = 0.3, Ao = 0.7, Qb = 0.02 
and h — 0.7 ^ has been shown to reproduce many of the 
properties of galaxies in the local Universe, such as their lu- 
minosity functions, the slope and scatter of the TuUy-Fisher 
relation, colours, sizes and metallicities (Cole et al. 2000) 
and also the clustering of galaxies in real and redshift space 
(Benson et al. 2000a,b). 

The model of Cole et al. (2000), like most other semi- 
analytic models of galaxy formation (e.g. Kauffmann et al. 

* We define Hubble's constant to be Ho = lOO/i kms~^Mpc~^. 
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1993; Somerville et al. 1999), includes a prescription for 
feedback due to energy input from supernovae and stellar 
winds. This is assumed to reheat cold gas and eject it from 
galaxies, thus inhibiting galaxy formation in low mass dark 
matter halos. Several studies of how this feedback may phys- 
ically operate can be found in the literature (Dekel & Silk 
1986; Mac Low & Ferrara 1999; Goodwin, Pearce & Thomas 
2000). This feedback is required in CDM models in order 
to produce a faint-end slope of the local galaxy luminosity 
function which is as shallow as that observed (White & Rees 
1978; Cole 1991; White & Frenk 1991), and also to produce 
galactic disks of sizes comparable to those observed (Cole et 
al. 2000). 

While ejection of gas by supernovae driven outflows is 
undoubtedly an important process (e.g. Martin 1999), other 
processes may also inhibit galaxy formation, for example, 
preheating of the IGM (Blanchard et al. 1992; Evrard & 
Henry 1991; Kaiser 1991; Valageas & Silk 1999), heat- 
ing of the gas inside galaxy and cluster halos (Wu, Fabian 
& Nulsen 2000; Bower et al. 2001), and the effects of a 
photoionizing background. The last of these is perhaps the 
best studied (see, for example, Efstathiou 1992; Thoul & 
Weinberg 1996; Katz, Weinberg & Hernquist 1996; Bullock, 
Kravtsov & Weinberg 2000). A photoionizing background 
both supplies heat to the gas through ionization, and reduces 
the rate at which the gas can cool by reducing the abundance 
of neutral atomic species which can be coUisionally excited. 
It thus raises the IGM temperature and so prevents it from 
collapsing into small halos, and also reduces the cooling rate 
of gas within halos and so reduces the fraction of baryons 
which can collapse to form a galaxy. For the formation of 
the very first objects at high redshift, cooling of the gas by 
molecular hydrogen is probably important, and one needs to 
consider the dissociation of these molecules by non-ionizing 
UV radiation (e.g. Ciardi et al. 2000), but these processes 
are only important well before the epoch of reionization in 
our model, since conversion of only a very tiny fraction of 
the baryons into stars is sufficient to produce enough UV 
radiation to dissociate all of the H2 molecules. 

In this section, we describe how we modify the model 
of Cole et al. (2000) to calculate the evolution of the IGM 
temperature and ionizing background, the suppression of gas 
collapse into halos by the IGM pressure, and the suppres- 
sion of cooling within halos by the ionizing background. Our 
modelling of the dynamical evolution of satellite galaxies 
within larger halos is described in 

2.1 Evolution of the Ionizing Background and the 
IGM Temperature 

We will treat the IGM as a mixture of six species (e. Hi, 
Hll, Hel, Hell and Hem) which interact with each other and 
with a uniform background of radiation emitted by stars 
and quasars. Since we are here primarily interested in the 
properties of low-redshift galaxies, we will not include H2 in 
our calculations, since it will be dissociated at high redshifts 
(e.g. Ciardi et al. 2000). We follow the evolution of the 
abundances of these species and the gas temperature for 
parcels of gas spanning a wide range in density contrast. The 
density contrast of eac h parcel is allowed to change with time 
as described in 



2.1.1 



Here, we treat all gas in the Universe 
as being part of the IGM. Since the fraction of the total gas 



content of the Universe which becomes part of a galaxy in 
our model is always small this is a reasonable approximation. 
Some gas should of course fall into the potential wells of dark 
matter halos (see §2.2). Since this gas typically occupies a 
small fraction of the volume of the Universe we ignore it for 
calculating the properties of the IGM. 

In the remainder of this section we describe in detail 
how we model the evolution of the ionizing background and 
IGM temperature. 



2.1.1 Evolution of Gas Density 

We wish to calculate the thermodynamic behaviour of gas 
in the IGM up until the point at which it falls into a viri- 
alised dark matter halo. The gas in the IGM will have a 
range of overdensities resulting from the growth of density 
fluctuations due to gravitational instability (we do not con- 
sider here the possibility of a multiphase medium which may 
also produce variations in gas density) . Since recombination 
rates, and consequently heating and cooling rates, depend 
on the gas density, it is necessary to take this evolving dis- 
tribution of densities into account in our model. 

We characterise the evolving distribution of gas den- 
sities via the probability distribution function (PDF), 
Pv(A,t), defined such that Pv(A,t)dA is the fraction of 
volume in the universe occupied by gas with a density con- 
trast A = p/p at time t, where p is the gas density at a point 
and p is the mean gas density in the Universe. Normalisation 
of this function to give the correct mean density and total 
mass requires that 



Pv(A,t)dA = 1, 



/o 
and 



APv(A,t)dA = 1. 



(1) 



(2) 



The fraction of mass with density contrast < A is given by 

(3) 



F{A,t) = f A'Pv(A',t)dA'. 
"'0 



We assume that as the gas density field evolves, the ranking 
of gas elements by density remains the same. The density 
contrast at time t of a gas element which has density contrast 
Ao at time to is therefore given by the solution of 



P(A[t],t) = P(Ao,to), 



(4) 



We can use eqn. (y) to calculate the evolution in over- 
density A[t] of individual parcels of IGM gas having differ- 
ent values of Ao, once the functional form and evolution of 
Pv(A,t) have been specified. In our standard model, we as- 
sume that the PDF has a log-normal form, which has been 
found to provide a reasonable description of the density dis- 
tribution produced by gravitational instability in the mildly 
non-linear regime (e.g. Coles & Jones 1991), 



Pv(A) 



exp 



(InA-lnA)^ 



2al 



(5) 



Here, ua determines the width of the distribution and the 
constants A and In A are fixed from the normalisation con- 
ditions (eqns. ^ and ^). The value of cta as a function of 
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time can be chosen to reproduce a desired baryonic clump- 
ing factor 



ydump 



A^Pv(A,t)dA, 



(6) 



where the overbar denotes a volume average. In particular, 
we will choose cta to reproduce the baryonic clumping fac- 
tor, /c™4'p"'^'^'' , derived by Benson et al. (2001a). In their 
calculation, Benson et al. (2001a) assumed that gas in the 
IGM essentially traces the dark matter except that pressure 
prevents the gas from falling into dark matter halos with 
virial temperatures less than lO'^K. They then calculated 
/ciump"'^'^' ~ 1 + cr^, where is the variance of the dark 
matter density field in spheres of radius equal to the radius 
of a 10*K halo (a^ was calculated from the smoothed non- 
linear dark matter power spectrum obtained using the tech- 
niques of Peacock & Dodds 1996). In the present work, the 
halo mass below which gas accretion is negligible varies as a 
function of time. Nevertheless, our estimate of the clumping 
factor should still provide a reasonable approximation. We 
note that at the redshift appropriate for Hi reionization in 
this work (see § [l.l| ) the two different clumping factors con- 
sidered by Benson et al. (2001a) are in fact very similar (see 
their Fig. 9). 

Once the evolution of the clumping factor has been cho- 
sen, our model results are insensitive to the particular func- 
tional form chosen for the PDF. For example, if instead of 
the lognormal distribution we use the form 



Pvi^) = (^) exp 



(|lnA-lnA|)^ 



2ai 



(7) 



which falls off much more rapidly away from In A = In A, 
this makes negligible difference to the evolution of the mean 
IGM temperature, ionization state and the spectrum of the 
ionizing background. We truncate the distribution of gas 
densities above A = 300, which is roughly the mean density 
contrast of halos at 2; = in our adopted cosmology, because 
reaction rates become extremely rapid for higher densities, 
making solution of the rate equations numerically difficult. 
Gas at higher overdensities accounts for only a small fraction 
of the total volume, and we have checked that moving the 
truncation point to larger A makes little difference to our 
results. 



2.1.2 Background Radiation 

We follow the proper number density of photons per unit 
frequency, n^, which evolves with time as 



dt 



a L 

^ ^ ^ ^ CO" u^i fv,j '^iij ; 



+ 5'. 



(8) 



where c is the speed of light, the term —3{d/a)ni, on the 
right-hand side represents the dilution of the number density 
by the Hubble expansion, and the term {a/a)d {un^) jdv de- 
scribes the effect of the redshifting of the photon frequencies. 
Here is the emissivity (i.e. number of photons emitted per 
unit volume, per unit time, per unit frequency), (t^,^ is the 
photoionization cross section for species i (Hi, Hel, Hell), 
/v,j is the fraction of the volume of the universe occupied 



by gas in density bin j, and n^j is the abundance of species 
i in density bin j. 

The photon number density is related to the background 
intensity by 



47r 



(9) 



where is the intensity per unit solid angle per unit fre- 
quency and hp is Planck's constant. 



2.1.3 Rate Equations 

The evolution of the abundances of the different ionization 
states of H and He is described by equations of the form 

= [ai{TiGM)ni+\nc — ai-\{TiGM)nine 

—Fc,i{TiGM)ninc + TcA~i{TiGM)ni^inc 

r^^jTlj -\~ T-f — 1 ^T-i— 1 ] 

1 dA(a) 



+ 



A(a) dt a 



3- I Ui 



(10) 



where for each atomic species H or He, i refers to the ioniza- 
tion state (i.e. i = 1, 2 for Hi, Hll and i = 3, 4, 5 for Hel, Hell, 
Hem), Tii is the proper number density, Tigm is the temper- 
ature, Qi is the recombination rate coefficient to i, Fc,! is the 
collisional ionization rate coefficient from i and F-^.i is the 
photoionization rate coefficient from i. The evolution of the 
electron density then follows from the conservation of the 
total number of electrons. 

We consider the evolution of a parcel of gas of density 
contrast A(t), which has a thermal energy per unit volume 
given by i5 = IkeTiGMn-tot, where ntot is the total number 
of particles per unit volume. The energy changes due to adi- 
abatic expansion/compression and atomic heating/cooling 
processes. Thus the evolution of E may be written as 



dE 
dt 



1 dA(a) 
A (a) dt 



E 



A' 



(11) 



where the first term represents adiabatic expansion or com- 
pression and the second represents atomic heating and cool- 
ing processes. E""" is the rate of heating per unit volume due 
to all heat sources (i.e. photoionization and Compton heat- 
ing) and is the rate of cooling per unit volume due to 
all heat sinks (i.e. Compton cooling and various atomic pro- 
cesses). We use the notation and A""^ to indicate rates 
of thermal energy gain/loss, as distinct from the usual ra- 
diative cooling function A which includes the entire energy 
of the photons emitted by recombinations. Note that the 
evolution of the gas density is entirely determined by the 
functional form of the PDF Pv(A) and the redshift evolu- 
tion of /clump and is unaffected by any heating/cooling of the 
gas. In reality, the gas density distribution should respond 
to differences in gas pressure. However, the effects of pres- 
sure forces should be important only on scales smaller than 
the Jeans length (or, more precisely, the filtering length to 
be introduced in the next subsection) , which always remains 
small ( <.lh~^Mpc) relative to the much larger scales over 
which we calculate volume averages. 

From eqn. (^) and the definition of E, we obtain the 
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following equation for the evolution of the IGM temperature: 



1 dTiGM 
TiGM At 



2 dA 
■'a 3A"dF 



-2- + — — + 3 



a^ntot At 



(12) 



The final term accounts for the effects of changes in the to- 
tal particle number density due to ionization/recombination. 
For a homogeneous IGM (A = 1) with no heating or cool- 
ing and no ionization or recombination, we have simply 



ritot oc a ^ and Tigm oc a ^. 



Equations |8|-|12| describe the evolution of a parcel of gas 
of specified final density contrast Aq. These equations, along 
with those describing the evolution of the background radia- 
tion spectrum, are solved for a range of Aq using a modified 
Bulirsch-Stoer method which is applicable to this stiff set of 
equations (Bader & Deuflhard 1983). The matrix decom- 
position that must be carried out as part of this method is 
efficiently achieved using a suitable sparse matrix package. 

The initial conditions for the abundance of each species 
and for the temperature are taken from the recfast code 
(Seager, Sasselov & Scott 2000) which accurately evolves 
the IGM through the recombination epoch (we typically be- 
gin our own calculations ai z = 200, at which point recom- 
bination is essentially complete but no significant sources of 
radiation have appeared in our model). 

We take photoionization cross sections from Verner et 
al. (1996), recombination rate coefficients from Verner & 
Per land (1996) and Arnaud & Rothenflug (1985) and col- 
lisional ionization rates from Voronov (1997). The cooling 
rate due to coUisional excitation of Hi was taken from Scholz 
& Walters (1991), while that for Hell was taken from Black 
(1981) with the modification introduced by Cen (1992) at 
high temperatures. The cooling rate due to free-free emission 
was computed using the Gaunt factors given by Sutherland 
(1998). 



2.2 Critical Mass for Collapse 

If the IGM has a non-zero temperature, then pressure forces 
will prevent gravitational collapse of the gas on small scales. 
In the absence of dark matter, the effects of pressure on the 
growth of density fiuctuations in the gas due to their self- 
gravity are described by a simple Jeans criterion, such that 
density fiuctuations on mass scales below the Jeans mass 
Mj are stable against collapse. However, this simple crite- 
rion needs to be modified in the case of non-linear collapse 
of the gas in the presence of a gravitationally dominant cold 
dark matter component which collapses to form dark matter 
halos. Gnedin (2000b) has obtained an analytical descrip- 
tion of the effects of gas pressure in this case, based on earlier 
work by Gnedin & Hui (1998). Using a linear perturbation 
analysis, Gnedin & Hui (1998) found that growth of density 
fluctuations in the gas is suppressed for comoving wavenum- 
bers fc > fcp, where the critical wavenumber fcp is related to 
the Jeans wavenumber kj by 



and fcj is defined as 



At'a\t' 



kj{t' 



dt" 



, a^{t") 



(13) 



^47rGpt, 



3/imH 
5kBT'IG^ 



1/2 



(14) 



In the above, ptot is the mean total mass density including 
dark matter, D{t) and H{t) are the linear growth factor and 
Hubble constant respectively as functions of cosmic time t, 
and ■ represents a derivative with respect to t. This expres- 
sion for fcp accounts for arbitrary thermal evolution of the 
IGM, through kj{t). Corresponding to the critical wavenum- 
ber fcp there is a critical mass Mp which we will hereafter 
call the filtering mass, defined as 



Mp = (47r/3)ptot(27ra/fcp)=' 



(15) 



The Jeans mass Mj is defined analogously in terms of fcj. 
In the absence of pressure in the IGM, a halo of mass A/tot 
would be expected to accrete a mass {flh/^o)Mtot in gas 
when it collapsed. Gnedin (2000b) found that in cosmologi- 
cal gas-dynamical simulations with a photoionized IGM, the 
average mass of gas Mgas which falls into halos of mass Mtot 
can be fit with the formula 



[1 + (21/3 _ l)MF/Mtot]3 ^^^^ 

with the same value of Mf as given by equations jl^ ) and 
(p^. The denominator in the above expression thus gives the 
factor by which the accreted gas mass is reduced because of 
the IGM pressure. Specifically, Mp gives the halo mass for 
which the amount of gas accreted is reduced by a factor 2 
compared to the case of no IGM pressure. 

In our model, we calculate the filtering mass Mp(2:) 
from equations (^^ and (p^), using for the IGM temperature 
the volume- averaged value Tigm ~ X/ fv.jTicM.j (where 
TiGM.i is the temperature of IGM gas in density bin j). 
We then modify the galaxy formation model of Cole et al. 
(2000) such that as each halo forms, it accretes a mass 
of hot gas given by equation (p^), multiplied by a factor 
(1 ~ /gal), where /gai is the fraction of the total baryonic 
mass (Qb/^^o)Mtot associated with that halo which has al- 
ready formed galaxies earlier on in progenitor halos. The gas 
which would have been accreted in the absence of IGM pres- 
sure is assumed to remain in the IGM, but is available for 
accretion later on in the merging process when another new 
halo is formed. The gas which does accrete is distributed 
within the halo as described by Cole et al. (2000). 



2.3 Cooling Rate of Gas in Halos 

The cooling of gas within dark halos, which controls how 
much of the gas can collapse to form galaxies, is also affected 
by the ionizing background. The gas within dark halos is at 
much higher densities than in the IGM, so we assume that 
it is in ionization equilibrium under the combined effects of 
atomic collisions and the external photoionizing background. 
We also assume that the halo gas is optically thin to the ion- 
izing background and to its own emission. While in our mod- 
els the mean metallicity of the IGM remains low enough that 
it has negligible effect on the cooling, this is not true for all of 
the gas in halos, some of which becomes significantly metal- 
enriched due to ejection of gas from galaxies by supernova 
feedback. We therefore use the photoionization code map- 
pings III, an updated version of the mappings ii code used 
by Sutherland & Dopita (1993), to calculate the radiative 
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cooling rate of gas in halos in collisional and photoioniza- 
tion equlibrium, including the effects of metals. Using this 
code, we calculate the net cooling/heating rate of the gas as 
a function of density, temperature, metallicity and redshift, 
using the photoionizing background predicted by our model 
at that redshift. We also include Compton cooling due to 
free electrons scattering off microwave background photons. 

Figure ^ shows the net cooling rate (i.e. the difference 
of heating and cooling rates) as a function of temperature, 
for gas in the presen ce o f the ionizing background from our 
fiducial model (see §4^), for a metallicity Z = 0.3Zq, at 
three different redshifts z = 0, 2 and 4. The cooling rates 
per unit volume are divided by n^, and calculated at den- 
sities nn = 1.3 X 10"^, 3.5 x 10"^ and 1.6 x lO^^cm"^ at 
redshifts 2 = 0,2 and 4 respectively, which correspond to the 
mean densities of gas in dark matter halos at those redshifts. 
We also plot the cooling curve in the absence of an ionizing 
background (dotted line). For z = and 2 = 2, gas cooler 
than T « 10^ '^K is actually heated rather than cooled in the 
presence of the ionizing background. We see that at 2 = 4, 
the cooling rate for gas at the average density for virialized 
halos is almost indistinguishable from the case of zero ioniz- 
ing background, a consequence of the high gas density at this 
redshift, and cooling is effective down to T « lO^K. On the 
other hand, at the lower redshifts plotted, photoionization 
almost completely suppresses cooling at T ^4 x 10*K. For 
gas at the halo virial temperature, the latter corresponds to 
a halo circular velocity of approximately 30km/s. 

In our model of galaxy formation, the gas in a dark 
matter halo is assumed to be isothermal at the virial tem- 
perature Tvir of the halo, and to have a uniform metallicity 
.^haio. The virial temperature is defined in terms of the cir- 
cular velocity 14 at the virial radius of the halo as 
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(17) 



At each timestep in our calculations we compute the age of 
the halo and the cooling time, defined as 



'^cool 



A(nH,Tvir, ^halo, z) 



(18) 



Equating TcooI to the age of the halo, we solve for the density 
of the gas which is just able to cool, and hence for the cooling 
radius, using the assumed density profile of the halo gas. 
We then calculate the mass and angular momentum of gas 
cooling in that timestep in the way described by Cole et al. 
(2000). Q 



2.4 Comparison with Numerical Simulations of 
the IGM 

Our model of the IGM is highly simplified, but we only re- 
quire it to predict a few volume averaged quantities, namely 
the IGM temperature and the spectrum of the ionizing back- 
ground. The advantage of our approach is one of speed, 
allowing rapid exploration of many different models. The 



^ Equation 4.3 of Cole et al. (2000) contains a typographical 
error — the factor /xmn should appear in the numerator, not the 
denominator. 
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Figure 1. The net cooling/heating function for gas at different 
redshifts in the presence of t he photoionizing background pre- 
dicted in our fiducial model (§1.2). We plot the absolute value of 
the cooling - heating rate per unit volume, divided by , for gas 
with metallicity Z = O.SZq, at redshifts 2 = (solid line), 2 = 2 
(dashed) and 2 = 4 (dot-dashed). At each redshift, we choose 
the gas density corresponding to the mean density in virialized 
halos at that redshift (thus nn = 1-3 X 10"'^, 3.5 X 10"-^ and 
1.6 X 10~^cm~'^ for 2 = 0, 2, 4 respectively). The dotted line 
indicates the cooling curve when no photoionizing background is 
present. The 2 = 4 curve is almost indistinguishable from this 
case. At low temperatures (T Ri lO'^ '^K) the 2 = and 2 = 2 
curves show a discontinuity, below which there is net heating 
rather than cooling 



disadvantages, compared to N-body /gas-dynamical simula- 
tions, are that it does not include the effects of spatial varia- 
tions in the ionizing background (no radiative transfer), and 
only includes the effects of gas density variations in a very 
approximate way. These limitations are likely to be most 
important just prior to full reionization, when there may be 
large spatial variations in the ionizing background and in 
the ionization state and temperature of the IGM. However, 
in this paper we are interested chiefly in calculating how an 
ionizing background suppresses galaxy formation, and these 
suppression effects only become strong after the IGM has 
been reionized, when our approximation of a uniform ioniz- 
ing background should be more accurate. 

We have tested the effects of the approximations in our 
model for the evolution of the IGM and ionizing background 
by comparing it to the N-body/gas-dynamical numerical 
simulations of Gnedin (2000a), which include the effects of 
the detailed spatial distribution of gas and ionizing sources, 
as well as an approximate treatment of radiative transfer. 
To do this test, we input into our model the same volume- 
averaged stellar emissivity and spectrum as measured from 
the simulations, assuming also the same cosmological pa- 
rameters. Fig. y compares predictions of our model with the 
same quantities measured from the simulations. Note that 
Gnedin's simulation was stopped at 2 = 5, so we cannot 
make any comparison at lower redshift. The left-hand panel 
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in Fig. ^ compares the volume- averaged IGM temperature, 
Hi and Hll fractions, and Jy(912A) as functions of redshift, 
and also the background radiation spectrum at z; = 9 (just 
prior to reionization for this model). The right-hand panel 
compares the Jeans and filtering masses predicted by our 
model with the values measured from the simulations by 
Gnedin (2000b). In the simulations, the filtering mass was 
determined by measuring the gas masses accreted by differ- 
ent halos and fitting these with the formula (p^. 

Overall, the level of agreement between the two ap- 
proaches is very good, although there are some differences 
in detail. The temperature of the IGM rises earlier in our 
model, and reionization occurs slightly earlier, presumably 
because the recombination rate in Gnedin's simulations is 
initially very high due to the ionizing sources forming in the 
highest density regions. The ionizing background at z ~ 9 
is in reasonable agreement with that from the simulation, 
although slightly higher. At wavelengths longwards of 912A 
our model predicts a significantly higher background. Here 
the gas is optically thin, so the details of absorption and the 
distribution of Hi are unimportant. It seems therefore that 
the approximate radiative transfer used by Gnedin (2000a) 
somewhat underestimates the background in the optically 
thin case. This is also apparent in the bottom right panel 
where we show J^(912A) as a function of redshift. Prior to 
reionization the two models predict very similar values, but 
afterwards our model reaches a significantly higher value 
than does Gnedin's. In the right hand panel of Fig. ^ we 
compare the Jean's and filtering masses. The Jean's mass in 
our model begins to increase sooner than in Gnedin's simu- 
lations (as expected from the earlier temperature rise in our 
model), and this difference is reflected in the filtering mass. 
Nevertheless, our simple model of the IGM reproduces with 
reasonable accuracy the evolution of the filtering mass in 
the numerical simulation. For our purposes, this is the most 
important result of the comparison, because the largest ef- 
fect of photoionization on galax y fo rmation is through the 
filtering mass, as we will see in §4.2. 



3.1 



MODEL FOR THE DYNAMICAL 
EVOLUTION OF SATELLITE GALAXIES 

Model for dynamical friction and tidal 
stripping 



When two dark matter halos merge, a new combined dark 
halo is formed. The largest of the galaxies they contained is 
assumed to become the central galaxy in the new combined 
halo, while the other galaxies become satellite galaxies in 
the new halo. These satellites evolve under the combined 
effects of dynamical friction, which makes their orbits sink 
towards the centre of the halo, and tidal stripping by the 
gravitational field of the host halo and central galaxy, both 
of the dark matter halos originally surrounding the satellites 
and of the stars they contain. The Cole et al. (2000) model 
included the effects of dynamical friction on the evolution of 
satellites, but did not include any treatment of tidal strip- 
ping. Since we are now interested in a more detailed study of 
the properties of the satellites around galaxies like the Milky 
Way (Benson et al. 2001b), we must improve our original 
model to include tidal effects on satellites. We do this by 



following the approach of Taylor & Babul (2000) (with a 
few modifications), following the orbits of satellites within 
host halos and making simple analytical estimates of tidal 
effects (both "static" tidal limitation and tidal "shocks"). 
Taylor & Babul (2000) show that this simple model for the 
evolution of satellite halos is able to reproduce well many of 
the results seen in high-resolution N-body simulations. We 
describe this part of our model briefiy, referring the reader 
to Taylor & Babul (2000) for a detailed discussion, but will 
highlight the differences between our model and theirs. 

We calculate the evolution of the orbit of each satellite 
galaxy in its host halo, under the influence of dynamical fric- 
tion and tidal stripping. We specify the initial energy E and 
angular momentum J of the orbit (after the halo merger) in 
terms of the parameters B!^ / -Rvir.host and e = J/ Jc respec- 
tively, where Rc{E) is the radius of a circular orbit with 
energy E, -Rvir.host is the virial radius of the host halo, and 
Jc (E) is the angular momentum of a circular orbit with en- 
ergy E. We choose a constant value of -Rl? /-Rvir.host for all 
satellites. Our standard choice is -Rc /^?vir,host = 0.5, which 
is representative of the median binding energy of satellite 
halos seen in high resolution N-body simulations (Ghigna et 
al. 1998) at the output time of the simulation. The typical 
value of -Rc /-Rvir.host for satellites just entering their host 
halo should presumably be somewhat higher, since by the 
output time satellites will have lost some energy through dy- 
namical friction. Lacking a direct measurement of the initial 
-Rc /-Rvir.host from simulations, we will simply use 0.5 as a de- 
fault, but will also explore other values to assess the impact 
of the uncertainty in this parameter on our final results. We 
select a value for the initial orbital circularity, e — J/ Jc by 
drawing a number at random in the range 0.1-1.0, which is 
a reasonable approximation to the distribution of circular- 
ities found by Ghigna et al. (1998). These choices for the 
initial orbital energy and angular momentum are the same 
as those of Bullock, Kravtsov & Weinberg (2000). Given the 
energy and angular momentum of the orbit, we determine 
the apocentric distance and begin integration of the orbit 
equations at that point, where tidal forces are weakest. 

We model the dark matter in both the host and satellite 
halos as an NFW density profile (Navarro, Frenk & White 
1997), modified by the gravity of the galaxy which has con- 
densed at the halo centre (the calculation of this adiabatic 
compression of the halo is described in detail by Cole et al. 
2000). The galaxy at the centre of each halo is modelled as 
a combination of disk and spheroid. The disk has a density 
distribution given by 



Pd{x,y,z) = pd.oexp 



sech' 



(19) 



where is the disk radial scale length, and h is the ratio of 
vertical to radial scale-length, which we take to be constant 
and equal to 0.1. The spheroid is modelled as a spherically 
symmetric r^''''-law. The masses and sizes of these compo- 
nents are determined as described by Cole et al. (2000). 

The satellite galaxy-|-halo moves under the infiuence of 
two forces. The first is just the net gravitational force due to 
the host halo and its central galaxy. The force due to the disk 
is calculated using the method of Kuijken & Gilmore (1989) . 
The second force is that due to dynamical friction, which 
we estimate using Chandrasekhar's formula (e.g. Binney & 
Tremaine 1987, section 7.1) 
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Figure 2. A comparison of properties of the IGM and tlie pliotoionizing background in our model (heavy lines) and in the numerical 
simulations of Gnedin (2000a) (thin lines). In our model we have assumed the same stellar emissivity as in the simulations. In the left- 
hand panel we show the volume averaged temperature (top left), the fractions of neutral (solid lines) and ionized (dashed lines) hydrogen 
(top right), the spectrum of the ionizing background at z ^ 9 (bottom left) and the evolution of Jy(912A) with redshift (bottom right). 
The right hand panel compares the Jean's and filtering masses from our model with the simulation. Dashed lines show the Jean's mass, 
while points with error bars show the filtering mass measured from the simulations and the solid line shows that predicted by our model. 



Fdf., = -47rG'Mf lnA.p,B(x)--^ 

|Vrcl,i 



(20) 



where InA^ is the Coulomb logarithm, pi the local density. 



B{x) = erf(a;) — 2xexp{—x )/\/Tr, x 



Vrcl 



\/V2o 



the velocity dispersion and Vrci,i is the relative velocity of 
the satellite and component i. We consider two components 
which contribute dynamical friction forces, namely the dark 
matter of the host halo and the spheroid of the central 
galaxy (which we treat together and indicate hereafter by 
a subscript "h"), and the disk of the host halo galaxy (indi- 
cated by a subscript "d" ) . We adopt the same values for the 
Coulomb logarithms as did Taylor & Babul (2000) (namely 
2.4 for the dark matter/spheroid and 0.5 for the disk), which 
fit the results of N-body simulations well. Taylor & Babul 
(2000) discuss in detail the possible choices for the Coulomb 
logarithms. While previous semi-analytic models have often 
used In Ah = InMh/Ms, we prefer to use the same value as 
Taylor & Babul (2000) for this present work. The dynami- 
cal friction force depends upon the mass Ms of the satellite. 
We include in this mass that part of the satellite galaxy and 
its dark halo which has not yet been stripped by tidal forces. 

For the disk velocity dispersion, we take = Vc/\/2, 
where Vc is the rotation speed of the disk (computed for a 
spherically averaged disk), as did Taylor & Babul (2000). 
This results in an unrealistically high velocity dispersion 
when applied to the Milky Way (where the observed 1-D 
velocity dispersion is approximately 30 to 40 kms~^). We 
prefer to use the Taylor & Babul (2000) value at present, 
but find that using a lower value (o"d = 0.2K:) has almost no 
effect on the results presented in this paper (e.g. the galaxy 
luminosity functions of Fig. ^ are hardly affected by this 
change). For the dark matter/spheroid system we find ah 



by integration of the Jean's equation (assuming an isotropic 
velocity dispersion) 



d(phgg) 
dr 



GMh(r) 



Ph(r), 



(21) 



where Mh is the total (i.e. dark plus baryonic) mass within 
radius r of the host halo. We assume that that dark matter 
follows the NFW profile for all radii outside of the virial ra- 
dius. Cole & Lacey (1996) show that the velocity dispersion 
calculated in this way is in reasonable agreement with that 
measured in N-body simulations. 

At each point in the orbit, we calculate the "static" 
tidal limitation radius of the satellite galaxy -|-halo, rt- This 
is the radius where the gravitational force of the satellite 
equals the sum of the tidal force from the host halo plus the 
pseudo-force due to the satellite's orbit. 



GMJrt) 



_ _d_ / GMhjR) 

" di? I i?2 



rt, 



(22) 



where R is the distance from the centre of the host halo, 
Ms(rt) is the total mass within radius rt of the satellite, 
and Lu is the instantaneous angular velocity of the satellite. 
Note that the factor of uj^ is strictly accurate only for cir- 
cular orbits. Here we follow Taylor & Babul (2000) and 
include this term for all orbits. For the purposes of this cal- 
culation and that of Uh, the mass of the host halo disk is 
spherically averaged (the assumption under which eqn. ^ 
was derived). In all other calculations of satellite dynamics 
we use the density distribution of eqn. ( |ig|) to describe the 
disk. Equation (M) is valid under the assumption that the 
satellite is much smaller than the host halo, which is true for 
all but a very small fraction of satellites in our calculations. 
Weinberg (1994) has argued that mass loss may oc- 
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cur at smaller radii than suggested by the above expression, 
due to heating of the satellite by gravitational shocking as 
it passes near the centre of the host halo. We adopt the ap- 
proach of Taylor & Babul (2000) to estimate the effect of 
this tidal shocking, and refer the reader to that paper for 
a complete description of the method. Briefly, during each 
fast shock (i.e. any shock for which the timescale is less than 
the internal orbital period of the satellite at its half-mass ra- 
dius), we calculate the rate of heating by tidal forces. The 
energy thereby deposited in the satellite causes the satellite 
to expand, pushing some material beyond the tidal radius, 
and so allowing more material to be removed by tidal forces. 
We define the effective tidal radius rt as the radius in the 
original satellite density profile beyond which material has 
been lost. When the satellite has been heated, this effec- 
tive tidal radius will therefore be less than rt as defined by 
eqn. ([22|). We remove matter from the satellite in spheri- 
cal shells outside of the effective tidal radius in the heated 
satellite. Note that we leave the density profile of material 
inside the effective tidal radius unchanged, so that the max- 
imum circular velocity in the satellite remains unchanged 
until the effective tidal radius is reduced below the position 
of the peak of the rotation curve (i.e. 2.16rs for a pure NFW 
dark matter rotation curve, but some other value when the 
baryonic contribution is included). 

Mass beyond the effective tidal radius of the satellite 
is removed gradually on the shorter of the angular orbital 
timescale, 2ti/u) (which becomes the orbital period for cir- 
cular orbits), and the radial infall timescale R/v-r.- Taylor & 
Babul (2000) only considered the angular orbital timescale, 
which results in very low mass loss rates for satellites on 
nearly radial orbits (as sometimes occur if the dynamical 
friction force is strong). It is then simple to calculate the 
mass remaining in the satellite (including both dark matter 
and baryonic components) and to use this to calculate the 
dynamical friction force exerted on the satellite. 

The orbit equations are integrated until one of three 
conditions is met; 

(i) The final redshift (i.e. the redshift at which we are 
studying the galaxy population) is reached. In this case we 
calculate the remaining mass and luminosity of the satellite 
galaxy after tidal limitation. 

(ii) The host halo merges to become part of a new halo. 
The satellite halo then becomes a satellite of the new halo 
and is assigned a new orbit in that halo. We begin inte- 
gration of the orbit equations again, but starting with the 
previous value of the effective rt for the satellite. 

(iii) The satellite merges with the central galaxy (which 
we assume happens when the orbital radius, R, first reaches 
-Rmerge, which wc take to be the sum of the half-mass radii 

of the host and satellite galaxies^. In this case we add to 
the central galaxy of the host halo the remaining mass of 
satellite galaxy at the time when the pericentre of its orbit 
first passed within iZmerge (even though some of this mass 
may have been stripped off since that time, its orbit will 
carry it into the central galaxy in any case), using the rules 



•t The model results are insensitive to the exact definition of 
merger time, as once R reaches such small radii it decreases very 
rapidly to zero. 




t/Gyr 

Figure 3. An example of the evolution a satellite galaxy orbit. 
The satellite enters the host halo at t = and merges with the 
central galaxy of that halo after approximately 2.5 Gyr. The solid 
line shows the orbital radius of the satellite as a function of time, 
showing that the orbit is decaying rapidly due to the efi'ects of 
dynamical friction. The dotted line shows the remaining mass of 
the satellite. Note that the mass does not begin to decrease until 
the first passage through pericentre, as before this the tidal forces 
felt by the satellite are not strong enough to remove any mass. 
See text for more details. 



described by Cole et al. (2000) (and possibly triggering a 
burst of star formation). 

Stellar mass stripped from satellite galaxies is added to 
a diffuse stellar component of the host halo, but is not con- 
sidered further in our models. Any cold gas stripped from the 
satellites is added to the hot gas reservoir of the host halo, 
and so may be able to cool again at a later time. Finally, 
as the satellite galaxy orbits in the host halo it continues to 
form stars, which causes some of the cold gas mass of the 
galaxy to be ejected according to the supernova feedback 
prescription of Cole et al. (2000). In the case of satellites 
with shallow potential wells, this can significantly alter the 
mass of the galaxy along its orbit. Therefore the mass of this 
reheated gas is removed from the satellite halo during the 
orbit. 

In Fig. ^ we show an example of a satellite orbit calcu- 
lated using the above model. The host halo has a mass of 
2x10^'^ Mq , a virial circular velocity of 440km / s and con- 
centration (defined here as the ratio of virial radius to NFW 
scale radius) of 5.9. The same three quantities for the satel- 
lite when it was still a separate halo are 3.5 x 10^^h~^ Mq, 
340km/s and 5.6, respectively. However, in this example, the 
satellite halo has already lost mass while being a satellite in a 
progenitor of the current host halo, which is why in this plot 
it starts from a mass of 8x 10^^ Mq. In this plot, the time 
t is measured from when the host halo formed, and the satel- 
lite orbit begins at the apocentre (approximately 250/i~^kpc 
from the centre of the host halo). The orbit decays rapidly 
due to the effects of dynamical friction, so that the satel- 
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lite makes three orbits before merging with the galaxy at 
the centre of the host halo at t « 2.5Gyr. The mass of the 
satellite is seen to decrease most rapidly when the satellite 
is close to pericentre. Note also that until just before the 
first passage through pericentre, the mass of the satellite is 
unchanging, as before this time tidal forces are simply not 
strong enough to strip any mass from the halo. 

We do not attempt here to model changes in the den- 
sity profile of the satellite galaxy -|-halo within the tidal ra- 
dius - the profile of the unstripped material is assumed to 
remain as it was before any stripping occurred. Nor do we 
account for any changes in the global properties of a galaxy 
which has lost mass to tidal forces (i.e. the galaxy keeps 
the same scale lengths, star formation timescale etc. as it 
had before any mass loss occurred). Numerical simulations 
of satellites undergoing tidal interactions (e.g. Mayer et al. 
2001) demonstrate that satellite mass profiles are affected 
by tidal interactions. Typically, they find a reduction in the 
amount of stellar mass within a given radius for radii within 
but comparable to the tidal radius. For low surface bright- 
ness galaxies (LSBs), this reduction can be up to a factor 
of around 2, but the effect is much weaker for high sur- 
face brightness galaxies (HSBs). At small radii, the stellar 
mass within a fixed radius is often increased by tidal interac- 
tions (through the production of a bar). The rotation curves 
of the galaxies are more seriously affected (presumably be- 
cause the spherically distributed dark matter is less strongly 
bound than the stellar disk), often being reduced by a factor 
2 for LSB galaxies (and somewhat less for HSB galaxies). 

Our calculation of merging times improves upon the 
simple estimates previously used in many semi-analytic 
models (which have often used results for satellites orbit- 
ing in isothermal halos, with no tidal stripping). Ifowever, 
we find that on average our approach predicts comparable 
merging timescales for satellite halos to the simpler treat- 
ment in Cole et al. (2000), although some fraction of satel- 
lites are predicted to have extremely long timescales, as they 
lose so much mass through tidal stripping that dynamical 
friction forces become extremely weak. 

3.2 Comparison with N-body simulations 

For our present purposes, we are most interested in whether 
our model reproduces the abundance of satellite halos, (or 
sub-halos) in a host halo typical of Milky Way-like galax- 
ies. We therefore compare the predictions of our model for 
the number of satellite halos with the results from high 
resolution, dark-matter-only N-body simulations of Milky 
Way-like halos in CDM models. Fig. ^ shows the compar- 
isons with simulations of ACDM (f2 — 0.3, A = 0.7) by 
Klypin et al. (1999) and of SCDM {Q = 1) by Moore et 
al. (1999a). In the ACDM simulation, the host halo mass 
is Mhaio = 1.1 X 10^^, and the minimum resolvable sub- halo 
mass is 3.3 x 10^h~^ Mq (corresponding to 20 particles). 
The corresponding quantities in the SCDM simulation are 
1.0 X 1O^^^"^M0 and 1.6 x IO^/i'^Mq. We make the com- 
parison in terms of the sub-halo velocity function A'^(> Vc), 
defined as the cumulative number of sub-halos per host halo 
with circular velocities greater than Vc, where Vc is defined 
as the peak circular velocity of the sub-halo. 

To compare our results to those of the N-body simu- 
lations, we run our semi-analytic model without baryons. 



Each dark matter halo then has a pure NFW profile. In the 
case of satellite halos, the NFW profile is truncated beyond 
a radius vt determined by the combined effects of static tidal 
limitation and tidal shocking. For an untruncated NFW pro- 
file, the circular velocity peaks at r,„ = 2.1 Bra, where Vs is 
the NFW scale radius. We assume that the density profile 
of satellite halos is unchanged within rt, so the peak Vc is 
VNFw(''m) if rt > rm, and VNFw(n) otherwise. We choose 
a host halo mass at z = equal to the value in the simula- 
tion, run 300 different realizations of the halo merger tree, 
and then take the mean A'^(> Vc) averaged over these re- 
alizations. In our model, sub-halos are only completely de- 
stroyed when they merge into the centre of the host halo. 
Tidal stripping reduces the mass of a halo, but is assumed 
never to destroy it completely. In our semi-analytic model, 
we can resolve much lower mass halos than can be resolved 
in the N-body simulations. Since there can be a wide range 
of sub-halo masses at a given value of the sub-halo circu- 
lar velocity Vc, it is essential to take into account this dif- 
ference in mass resolution when we compare to the simu- 
lations. Therefore, to calculate A^(> Vc), we discard from 
the semi-analytic model all sub-halos with masses (within 
rt) smaller than the minimum resolvable sub- halo mass in 
the simulation. As Fig. ^ shows (compare the solid and dot- 
dashed curves), this mass cut produces a large reduction in 
the number of satellites below Vc = 20 — 40kms~^, and is 
very important for matching the simulation results. 

The upper left and right panels of Fig. ^ show the com- 
parison of our model with the ACDM and SCDM simu- 
lations respectively. The solid curves show the prediction 
for i?c/^vir,host = 0.5, iucludiug the mass resolution cut. 
The dotted lines on either side of the solid line show the 
10%-90% range of the distribution seen among the differ- 
ent realizations. This range is larger for the ACDM than 
SCDM model, which mainly results from the smaller num- 
ber of sub-halos per host halo in the former case. In the 
same panels, the upper and lower dashed lines show the 
effect of changing the assumed initial orbital energy to 
-Rc/i?vir,host = 0.75 and 0.25 respectively. We see that our 
standard value J?c/^vir,host = 0.5 gives significantly better 
agreement with the N-body simulations for both ACDM and 
SCDM. 

The lower left-hand panel of Fig. ^ shows in more de- 
tail the separate effects of dynamical friction, static tidal 
limitation (i.e. rt as defined by eqn. and tidal shocks 
on the velocity distribution, for the ACDM case. All of the 
curves plotted there assume i?c/-Rvir,host = 0.5, and include 
the cut in sub-halo mass corresponding to the resolution 
of the N-body simulation. The solid line includes all of the 
above processes, and is therefore identical to the solid line 
in the upper left-hand panel. The dotted line contains none 
of these processes, so sub-halos never merge and are never 
tidally stripped. Switching on dynamical friction results in 
the dashed line, which greatly reduces the number of high 
Vc (relatively massive) subhalos, but is much less important 
for the low Vc halos. Switching on static tidal limitation 
(and keeping dynamical friction switched on) results in the 
dot-dashed line. This greatly reduces the number of low Vc 
halos, as these are strongly affected by tidal forces, once 
the mass cut is included. The number of high 14 halos ac- 
tually increases somewhat, since tidal limitation is able to 
reduce the mass of these halos and so reduce the strength 
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Figure 4. The number of satellite lialos as a function of circular velocity from the semi-analytic model compared to N-body simulations, 
for Milky Way-like halos in CDM models. N{> Vc) is the cumulative number of sub-halos per host halo, with Vc defined as the peak 
circular velocity of the sub-halo. The simulations are of a ACDM model from Klypin et al. (1999) and of a SCDM model from Moore 
et al. (1999a). In each panel, the solid points with error bars show the N-body simulation results, while the lines show the semi- 
analytic predictions for different assumptions. Upper two panels: the light solid lines show the semi-analytic prediction (averaged over 
300 realizations), including in the semi-analytic model the same cut on sub-halo mass as in the simulations. The dotted lines on either 
side of the solid line show the 10%-90% range of the distribution around the mean value. The dot-dashed lines show the result from the 
semi-analytic model when no mass cut is applied. These results are for Rc/^ir,host = 0-5. The dashed lines show the results if instead 
we assume iic/^vir,host = 0-75 or 0.25, including tidal stripping and the mass cut. Lower panel: this shows the contribution of different 
physical processes in the semi-analytic model, for the case of ACDM. The dotted line shows the predicted sub-halo velocity distribution 
for the case of no dynamical friction, no static tidal limitation and no tidal shocking. Switching on dynamical friction produces the 
dashed line. Adding in static tidal limitation gives the dot-dashed line, and finally switching on tidal shocking produces the solid line. 
In this panel, all the curves are for R^/ Rvii:,host = 0-5, and all include the same cut on sub-halo mass as in the N-body simulation. 



© 0000 RAS, MNRAS 000, 000-000 



12 A. J. Benson, C. G. Lacey, C. M. Baugh, S. Cole,& C. S. Frenk 



of the dynamical friction forces which they experience. The 
remaining difference between the dot-dashed and solid lines 
is accounted for by tidal shocking. The overall effect of tidal 
stripping is to reduce the number of halos at low Vc by a 
factor ~ 10. 

An important difference between our model and the 
similar calculation by Bullock, Kravtsov & Weinberg (2000) 
is we can resolve sub-halos within sub-halos (i.e. we record 
all branches of the merger tree and so halos merging with 
the final halo may have substructure of their own) , whereas 
Bullock et al. considered only sub-halos (i.e. merging halos 
were assumed to have no substructure of their own). This 
distinction is important, since, in our model of galaxy for- 
mation, every branch of the merger tree can potentially host 
a galaxy (providing its virial temperature exceeds 10'*K and 
so is able to cool efficiently) . Hence, when no tidal stripping 
is applied we find many more satellites at a given Vc than 
did Bullock et al. However, when tidal stripping is included, 
both our model and that of Bullock et al. are in reasonable 
agreement with the N-body results (although Bullock et al. 
compared only with the ACDM simulations of Klypin et al. 
(1999)). 

We remind the reader that our choice of Rc/Rvir = 0.5 
was originally motivated by the measurement of the orbital 
energy distribution of all the satellites existing in a halo at 
the final output time of an N-body simulation. However, in 
our model, we use this R^/Rvir as the initial value for each 
satellite after it joins the main halo. Plausibly we should 
use a higher value, since when satellites first fall into a host 
halo they should be less bound than at any subsequent time. 
This would reduce the effectiveness of tidal limitation in 
our model (e.g. compare the curves for R°/R^ir = 0.50 and 
0.75 in Fig. ^). Also, as noted above, we do not include any 
adjustment in the density profile of the material within the 
tidal radius in response to stripping of material from larger 
radii. This would be expected to make satellites less bound 
and to enhance the process of tidal stripping, and also to 
lower the circular velocity. These two effects work in opposite 
directions, but it is not clear which is the dominant process. 
N-body suggestions suggest that the effect of the latter on 
the sub- halo peak circular velocities Vc is in fact fairly small; 
Ghigna et al. (2000) find in their high-resolution simulations 
that for sub-halos where the tidal radius is larger than the 
initial peak-Vc radius, Vc typically changes by only ~ 20% 
due to tidal effects. For now, we simply note that Rc/Rvir = 
0.50 does produce reasonable agreement with the numerical 
results, and so we adopt this throughout the remainder of 
this paper. 

Our semi-analytic model includes the effects of baryonic 
collapse on the mass profiles of the host and satellite halos, 
although this effect is turned off when we compare to pure 
dark matter N-body simulations. While baryonic dissipa- 
tion makes satellite halos more strongly bound, and so more 
resistant to tidal limitation, it also makes the tidal forces 
of the host halo and central galaxy stronger. The central 
galaxy of the host halo also contributes to satellite destruc- 
tion through its contribution to dynamical friction. In halos 
of mass ~ 10^^ Mq, the net result of including the bary- 
onic components is to further reduce the number of satel- 
lite halos (compared to a pure dark matter calculation). We 
find, for these halos and with our standard galaxy formation 
model (see §[l.2|) , that the number of satellites at a given Vc 



is reduced by around 40% at K 
60% at Vc = 15km/s. 



4 RESULTS 



60km/s, and by about 



We are now able to explore in a self-consistent way the effects 
of a photoionizing background on the properties of galaxies 
and also the effects of galaxies on the IGM. We will begin in 
§4.1 by obtaining a self-consistent model and explor ing the 
evolution of the IGM and ionizing background. In §4.2 we 
examine the effects on the global properties of galaxies in the 
fiducial model of Cole et al. (2000). The properties of satel- 
lite galaxies will be explored in a separate paper (Benson et 
al. 2001b). 



4.1 Properties of the IGM and Ionizing 
Background 

4.1.1 Star Formation History 

Our starting point is the fiducial model of Cole et al. (2000) , 
modified in the ways described in §^ and We use this to 
predict the star formation history and associated emissivity 
in ionizing photons as a function of redshift. We resolve all 
halos that are able to cool in the redshift interval to 25 
to ensure that all ionizing photons are accounted for. To 
determine the spectrum of emission from these stars, we 
tabulate the mean star formation rate per unit volume from 
our model as a function of both cosmic time and metallicity, 
d? pi{t,Z)/AtAZ. The stellar emissivity per unit volume at 
cosmic time t is then simply 



dt'dZ 



^{t~t',Z) dZdt', 



(23) 



where <E>(f, Z) is the spectral energy distribution of a stellar 
population of age t and metallicity Z, which we take from 
the models of Bruzual & Chariot (1999). 

To account for the effects of absorption by dust and 
gas in galaxies on the ionizing emissivity, we multiply the 
above expression (^sj) at wavelengths A < 9 12 A by a con- 
stant factor /osc, defined as the fraction of ionizing pho- 
tons produced by stars that escape through the dust and 
gas of the galaxy's interstellar medium (ISM). To calculate 
the effects of dust absorption for the non-ionizing radiation 
at A > 912A, we use the same approach as in Cole et al. 

(2000) . The value of /esc for the ionizing photons is uncer- 
tain, so we will present results for two values, /osc = 10% 
and /esc = 100%, which we believe to bracket a reason- 
able range. The value /esc ~ 100% results in our model in 
an emissivity in Lyc (A < 912A) photons from galaxies at 
z — 3 that agrees with the recent observational estimate 
for Lyman-break galaxies by Steidel, Pettini & Adelberger 

(2001) (after allowing for the differences in the assumed cos- 
mological models), and predicts reionization of hydrogen at 
2 ~ 8, compatable with measurements of the Gunn-Peterson 
effect in quasars. On the other hand, the value /esc = 10% 
is more consistent with both observational (e.g. Leitherer et 
al. 1995; Steidel, Pettini & Adelberger 2001) and theoreti- 
cal (e.g. Dove, ShuU & Ferrara 2000; Benson et al. 2001a) 
estimates of the escape fraction at both low and high red- 
shift. The observational estimate of L^(900A)/L^(1500A) 
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for the Lyman-break galaxies by Steidel, Pettini & Adel- 
berger (2001) implies /esc ~ 10 — 40%, allowing for uncer- 
tainties in dust extinction and in the emission at A < 912A 
predicted by stellar models. The reason why we require a 
larger /eac than Steidel, Pettini & Adelberger (2001) to pro- 
duce the same net ionizing emissivity at 2 = 3 is that the 
Cole et al. (2000) model predicts too low a typical 1500A 
luminosity for Lyman-break galaxies, once we include dust. 
Our model with fasc = 10% predicts an ionizing background 
that is in better agreement with observational estimates at 
z < 4.5, but also predicts a low redshift for reionization, 
z ~ 5, that is barely compatable with the Gunn-Peterson 
constraints. These issues are discussed in more detail below. 
Since it is the redshift of hydrogen reionization and reheat- 
ing of the IGM that appears to be the most important in 
determining the effects on galaxy formation, we will take 
/esc ~ 100% as our standard case. 

We also include the contribution to the ionizing emis- 
sivity from quasars, according to the observational parame- 
terization of Madau, Haardt & Rees (1999). Their parame- 
terization is based on fitting observational data on numbers, 
magnitudes and redshifts of quasars at z < 4.5, assuming an 
Einstein-de Sitter cosmology. To obtain the emissivity in our 
chosen cosmology, we must allow for the dependence of the 
observationally-inferred luminosities and number densities 
on the assumed cosmological model. We therefore use 



stellar 

Qiiasar 

Heavy lines = HI 
Thin lines = Hell 




Figure 5. The emissivities in Hi and Hen ionizing photons (heavy 
and thin lines respectively) per comoving volume as a function 
of redshift. Solid lines show the emissivity from stars assuming 
/esc = 100%, while dashed lines show that from quasars. 



e{z) = eMHR(z) 



dL(z) 
dr'(.) 



fdV{z}/dz^ 
\ dV{z)/dz 



(24) 



where eMHR(z) is the emissivity from Madau, Haardt & Rees 
(1999) (measured from their Fig. 2), d-L{z) is the luminosity 
distance and dV(z)/dz is the comoving volume per unit red- 
shift. Functions with superscript (EdS) are calculated in the 
Einstein-de Sitter cosmology, those without in the cosmol- 
ogy of our fiducial model. We use the same expression to ex- 
trapolate the quasar contribution to z > 4.5. We note that 
even at z < 4.5, the use of eMiiR{z) from Madau, Haardt 
& Rees (1999) involves a considerable extrapolation of the 
quasar luminosity function down to luminosities not directly 
observed. 

We then use the total (i.e. stellar plus quasar) ioniz- 
ing emissivity in calculating the thermal evolution of the 
IGM (and hence the filtering mass) and the ionizing back- 
ground, both as functions of redshift. The stellar emissivity 
must be determined self-consistently with the feedback ef- 
fects on galaxy formation from th e IG M pr essu re and ion- 
izing background, as described in 



2.3 



We do this 



jj2.2| and 

by means of an iterative procedure, starting from a galaxy 
formation model computed ignoring these feedback effects, 
calculating the ionizing background in this model, using this 
as input in calculating a revised model including the pho- 
toionization feedback effects, and repeating this cycle until 
we have a model whose star formation history is consistent 
with the photoionizing background that it produces. Cole 
et al. (2000) chose the parameters of their fiducial model 
to match certain observations of the local galaxy popula- 
tion, in particular the luminosity function of galaxies in the 
B and K-bands. We find that if we keep the same param- 
eter values as used by Cole et al. (2000), then when we 
include the photoionization feedback, our model still pro- 
duces an acceptable fit to these luminosity functions. The 



only change is a small adjustment of T (which determines 
mass-to-light ratios) from 1.38 in Cole et al. (2000) to 1.32 
(we adjust the recycled fraction in our chemical evolution 
model accordingly). As will be discussed in more detail in 
§4.2.1, the faint-end slopes are now somewhat flatter than 
before, but this is consistent with recent determinations of 
the luminosity functions. 

The net emissivity in Hi and Hell ionizing photons 
from stars and quasars for our standard model is shown 
in Fig. ^. For /esc = 100% stars always dominate the pro- 
duction of Hi ionizing photons, but quasars dominate the 
production of Hell ionizing photons until 2: « 1 when the 
rapidly falling quasar emissivity leads to stars becoming the 
dominant source. 

Figure ^ shows the star formation rate per unit comov- 
ing volume in our standard model (/esc = 100%) as a func- 
tion of redshift. The comparison of the solid line (final itera- 
tion) and dashed line (penultimate iteration) shows that the 
model has converged to a self-consistent star formation his- 
tory in the presence of the photoionization feedback.]^ The 
dotted line shows for comparison the star formation history 
from the model of Cole et al. (2000), with no photoion- 
ization feedback and scaled to the value of T used in our 
standard model. For 2 J>10 this is identical to that of our 
new model, a fact which is not surprising, since we use the 
same parameters as did Cole et al. (2000), and at these 
redshifts photoionization has yet to have much effect on the 
IGM. (Recall that gas in halos with Tvir ;^10'*K is assumed 



3 We restart our calculation of the star formation rate at several 
intervals in redshift to ensure all halos are resolved. Since the 
merger trees used in our model do not reproduce exactly the 
Press-Shechter mass function at an earlier redshift this results in 
small discontinuities in the star formation rate visible in Fig. ^. 
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Penultimate Iteration 
Cole et al, (2000) 
No SNe feedback 
No SNe feedback 
or photoionization 




Figure 6. The star formation rate per comoving volume in our 
standard model as a function of redshift. The solid line represents 
the star formation rate used to compute the temperature of the 
IGM and the evolution of the ionizing background. The dashed 
line indicates the star formation rate in the penultimate iteration 
of the model, indicating that convergence has been reached over 
the range of interest. The dotted curve indicates the star forma- 
tion rate in the model of Cole et al. (2000) scaled to the value 
of T in our standard model for the purposes of this comparison. 
We also show a model with the effects of supernovae feedback 
switched off, both with (dot-dashed line) and without (short- 
dashed-long-dashed line) the effects of photoionization feedback 
included. The epochs of Hi and Hen reionization in the standard 
model are marked by vertical dotted lines. The small discontinu- 
ities in the star formation rate arise as we recompute our model 
at several intervals in redshift to ensure all halos are resolved. 



to be unable to cool, even in the absence of an ionizing back- 
ground.) Beginning just before Hi reionization (as the IGM 
is being reheated) the star formation rate in our model falls 
below that of Cole et al. (2000) as the filtering mass rises. 
By z ~ 4 the star formation rate has recovered to the Cole 
et al. (2000) value as the continued formation of structure 
has created many halos well above the filtering mass, and it 
is these which contribute most to the star formation rate. 
(The filtering mass is growing only rather slowly during this 
period.) The reionization of Hell leads to a second episode of 
reheating, leading to an increase in the filtering mass which 
again suppresses star formation rates below the Cole et al. 
(2000) values. Once again, by z = the differences have be- 
come very small, as star formation becomes dominated by 
galaxies in halos well above the filtering mass. The effect 
is rather small however, with star formation rates being re- 
duced by around 25% at most. The reason why the effects of 
photoionization feedback on the star formation history are 
quite modest in our model is that we also include supernova 
feedback according to the prescription of Cole et al. (2000) . 
This greatly suppresses star formation in halos with circular 
velocities Vc ^ 200kms^^, which includes the range of halo 
masses that are also affected by photoionization feedback. 



For comparison, we have also computed a model in 
which the feedback from supernovae is completely turned off. 
The other parameters in this model are identical to those in 
Cole et al. (2000), apart from T, which is reduced to 0.95 to 
match the bri ght e nd of the present-day galaxy luminosity 
function (see §4.2.1). The star formation rate as a function of 
redshift in this model, with /esc = 100% and photoionization 
feedback turned on, is shown by the dot-dashed line in Fig. ^ 
The absence of feedback from supernovae results in a much 
higher star formation rate than in our standard model at all 
redshifts, but especially at high redshift, where star forma- 
tion is mostly occuring in small halos which are the most 
strongly affected by supernova feedback. As a result, reion- 
ization occurs significantly earlier in this model, ai z — 11.5 
for Hi and at z — 4.5 for Hell. The short-dashed-long-dashed 
line shows the star formation rate when feedback from su- 
pernovae and from photoionization are both turned off. The 
reduction in the star formation rate after reionization due 
to photoionization feedback is seen to be much larger when 
there is no feedback from supernovae, than in our standard 
model which includes supernova feedback. 



4-1.2 Thermal and Ionization History of the IGM 

Figure ^ shows various properties of the IGM in our standard 
model (with /esc = 100%). We also show selected results for 
models with a uniform IGM (i.e. /dump = 1) and for a model 
with a lower escape fraction, /esc = 10%. 

The top left-hand panel shows the volume averaged 
temperature of the IGM as a function of redshift. At around 
z = 20 ionizing photons from stars begin to heat the IGM 
(at higher redshifts the temperature scales as expected from 
adiabatic expansion, Tigm oc {1 + z)^). This results in the 
gas reaching a temperature of approximately lO^K at z ~ 10 
(somewhat before reionization of Hi), at which point atomic 
cooling processes balance the photoheating. The tempera- 
ture then decreases until z ~ 7, when the photoionization 
of Hen by emission from quasars leads to a second period 
of heating which lasts until 2 « 5. After this, the gas cools 
rapidly until z = 0, the cooling being due to adiabatic ex- 
pansion (due both to Hubble expansion and the expansion 
of the gas in voids). The redshifts of reionization of Hi and 
Hell in the standard model are marked by vertical dotted 
lines. (We define the redshift of reionization somewhat arbi- 
trarily as the point where 99% of the species in question has 
been ionized. Since reionization takes place rapidly the ex- 
act definition is not important.) With a uniform IGM, the 
temperature evolution is unchanged up to just before the 
reionization of Hi, since the clumping in the standard model 
is relatively small at these high redshifts and the volume- 
weighted mean temperature shown here is dominated by the 
contribution from gas close to the mean density. At lower 
redshifts the clumpy IGM cools more rapidly as adiabatic 
expansion of gas in voids cools the gas (and these regions 
are strongly weighted in the volume-averaged temperature) . 
The low /esc model heats the Universe later, as expected. 
There is little difference in the peak temperature reached, 
which is essentially fixed by atomic physics, and the late 
time temperatures are very similar to those of the standard 
model. 

The top-right hand panel shows the mean ionization 
state of hydrogen and helium in the IGM as a function of 
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Figure 7. Upper left-hand panel: The volume averaged temperature of the IGM as a function of redshift. The redshifts of reionization 
for Hi and Hell in our standard model are indicated by vertical dotted lines. The solid line shows the standard model (/esc = 100%), with 
the other lines showing the results for /dump = 1 and /esc = 10% as indicated in the panel. Top right-hand panel: The volume averaged 
ionization state of hydrogen and helium as a function of redshift. The quantity shown is ni/utot where Jitot is the total abundance of 
the element in question in all ionization states. The epochs of Hi and Hell reionization in our standard model are indicated by vertical 
dotted lines (we define these as the time at which Ui/utot reaches 0.99). We also show nni/itot for models with /dump = 1 and /esc = 1 
as indicated in the panel. Lower left-hand panel: The .leans mass (dashed line) and filtering mass (line types as defined in the panel) 
as a function of redshift. The filtering mass is shown for the standard model and also for models with /dump = 1 and /osc = 10%. The 
redshifts of reionization for Hi and Hell in our standard model are indicated by vertical dotted lines. Lower right-hand panel: The solid 
line shows the halo circular velocity (left-hand axis) or halo virial temperature (right-hand axis) that correspond to the filtering mass in 
the standard model. The dashed line shows the mean IGM temperature (repeated from the top left panel). 



redshift. The quantity plotted is the average fraction of hy- 
drogen or helium in each ionization state (e.g. for the HI 
fraction xm, we plot (nHii)/([nHi +nHii]), where () denotes a 
volume average) . Reionization is a much more rapid process 
than reheating (as has been noted previously by Gncdin & 
Ostriker (1997) and Valageas & Silk (1999) for example). 
Hei is ionized by stellar photons almost simultaneously with 



Hi at a redshift of 8, but Hen is not reionized until much 
later {z w 4) when the harder ionizing photons from quasars 
become abundant. Note that the initial decline in xm is sim- 
ilar in clumpy and uniform IGMs. As noted above, at these 
redshifts volurnc-avcragod quantities in the clumpy case are 
dominated by gas close to the mean density, so we do not 
expect much difference from the uniform case. Once started 
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though, reionization is completed much more rapidly in the 
case of a uniform IGM. In the clumpy IGM the comple- 
tion of reionization is delayed by the small fraction of high 
density gas, which is reionized last. With a low escape frac- 
tion of 10% reionization does not occur until much later, at 
z ~ 5.5 (which may be slightly too low to be consistent with 
recent measurements of the Gunn-Peterson effect in quasars 
at z ~ 6: Fan et al. 2000; Djorgovski et al. 2001; Becker et 
al. 2001), but otherwise proceeds in much the same way. 

The lower left-hand panel shows the evolution of the 
Jean's and filtering masses with redshift. Note that both of 
these are defined as halo masses, not baryonic masses. The 
Jean's mass simply tracks the temperature of the IGM, while 
the filtering mass approximately tracks the Jean's mass, but 
with a significant delay. As a result the filtering mass can 
be up to 1000 times lower than the Jean's mass during the 
first episode of reheating. However, at lower redshifts the 
two are much more comparable and by 2 = the filtering 
mass is around 60% of the Jean's mass. Note that the pe- 
riod of cooling from z = 9 to z — 7 (during which time Hi 
reionization has finished, but Hen reionization has yet to 
begin) the Jean's mass decreases slightly with time, and as 
a consequence the filtering mass grows only slowly. Using a 
uniform rather than a clumpy IGM affects Mp only at low 
redshifts, where the lack of cool void gas in the uniform case 
results in a slightly larger filtering mass. Although the filter- 
ing mass does not begin to rise until later in a model with 
/esc = 10%, it actually rises above the /esc ~ 100% model 
at late times since the IGM has actually been hotter in the 
recent past in this model. 

The lower right-hand panel shows the values of the halo 
circular velocity at the virial radius, and corresponding halo 
virial temperature (eqn.^^, that correspond to the halo fil- 
tering mass. Also shown for comparison is the average IGM 
temperature as a function of redshift. It can be seen that, ac- 
cording to the filtering mass prescription of Gnedin (2000b) 
that we use, the critical halo virial temperature below which 
baryonic collapse into halos is suppressed by 50% in mass 
can be much greater than the IGM temperature (by a factor 
60 at 2 = in our standard model). In our standard model, 
this temperature peaks at Tvir ~ lO^K, corresponding to 
Vc « 60kms~^, even though the IGM temperature is never 
significantly above lO^K. Clearly, it will be very important 
to test the accuracy of Gnedin's filtering mass prescription 
in greater detail using future high resolution simulations. 
However, we note that similar results for the halo circular 
velocity below which baryonic collapse is 50% suppressed 
were also found by Quinn et al. (1996), from SPH simu- 
lations, and Thoul & Weinberg (1996), using a ID hydro 
code (they both found Vc ~ 50kms~^ at z ~ 2). 



4.1.3 The Ionizing Background 

Figure |^ shows the evolution with redshift of the predicted 
ionizing background at the Lyman limit, J^{912K). This is 
compared to observational estimates from the proximity ef- 
fect in quasar spectra (e.g. Scott et al. 2000, and references 
therein), and upper limits from observational searches for 
Ha fluorescence from extragalactic Hi clouds at low redshift 
(Vogel et al. 1995) and for Lya fluorescence from Ly-limit 
clouds at high redshift (Bunker, Marleau & Graham 1998) . 
It can be seen that the background predicted for /esc = 10% 




Figure 8. Ionizing background vs redshift. The solid and dot- 
dashed lines are the predicted background intensity at the Ly- 
man limit, J„(912A), for /esc = 100% and 10% respectively. The 
dashed and dotted lines show the separate contributions from 
stars and quasars in the standard model with /esc = 100%. The 
rectangular boxes and the datapoint at z = 4.5 are observational 
estimates based on the proximity effect (Kulkarni & Fall 1993; 
Batjlik, Duncan & Ostriker 1988; Scott et al. 2000; Cooke, Es- 
pey & Carswell 1997; Williger et al. 1994). The upper limits are 
based on searches for Ho (Vogel et al. 1995, square; Weymann 
et al. 2001, star) or Lya (Bunker, Marleau &; Graham 1998) 
fluorescence; the latter limit is somewhat model-dependent. 



is reasonably consistent with observational estimates, while 
that predicted for /esc = 100% is 5-10 times too high. While 
this comparison seems to favour the model with /esc = 10%, 
the estimated ionizing emissivity of Lyman-break galaxies 
at 2 ~ 3 is only reproduced in our model with /osc ~ 100%. 
Furthermore, with /osc = 10%, reionization seems to oc- 
cur too late compared to observational constraints from the 
Gunn-Peterson effect, as already mentioned. This contradic- 
tion in part arises because our simple IGM model appears to 
underestimate the opacity of the IGM to ionizing photons at 
epochs when the IGM has been almost completely reionized. 
According to Madau, Haardt & Rees (1999), at 2 ^5, the 
Lyc opacity is dominated by the discrete absorbing clouds 
with neutral hydrogen column densities Nhi ~ lO^^cm"^, 
which produce the Lyman-limit absorption features seen in 
quasar spectra. Madau et al. estimate the opacity as a func- 
tion of redshift based on the observed statistics of quasar 
absorption lines, and find that the universe becomes opti- 
cally thin to Lyc photons only at z ^1.6. In contrast, our 
IGM model, which lacks these absorbing clouds, already be- 
comes optically thin to ionizing photons at 2 ~ 6 for the 
case /esc = 100%. This explains why our standard model 
produces a much larger ionizing background at 2 « 3 than 
Steidel, Pettini & Adelberger (2001) calculate from combin- 
ing their estimate of the ionizing emissivity of Lyman-break 
galaxies (which our model matches) with Madau et al.'s es- 
timate of the Lyc opacity. 

We would need to develop a much more sophisticated 
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IGM model in order to include the effect of discrete clouds 
on the Lyc opacity in a way that was both self-consistent 
and agreed with observations of quasar absorption lines. Ac- 
cording to our models, photoionization affects galaxy forma- 
tion primarily through the effect of the IGM pressure (which 
mainly depends on the redshift of reionization) rather than 
on the cooling within halos (which depends on the ionizing 
background at the redshift when the halo forms). Therefore, 
we believe that the defficiencies of our model as regards pre- 
dicting the ionizing background after reionization should not 
seriously affect the predictions that we make for galaxy for- 
mation. 



4.2 Effects on Galaxy Properties 

We now use our model, together with the properties of the 
IGM and ionizing background calculated in the previous 
subsection, to investigate the effects of photoionization feed- 
back on the global properties of galaxies at 2 = 0, and to 
briefly consider the effects on galaxies at higher redshifts. 



J^.2.1 Luminosity Functions 

In Fig. ^ we present the B and K-band luminosity functions 
for this model at z=0, and compare them to a selection of 
observational data. Note that all model galaxy luminosities 
include extinction by dust, calculated using the model of 
Ferrara et al. (1999) as described by Cole et al. (2000). 
The heavy solid line shows the result from our standard 
model, while the thin solid line shows that from the model 
of Cole et al. (2000) . Brighter than the galaxy luminosity 
function is mostly unaffected by the inclusion of the effects 
of photoionization. Fainter than this differences become ap- 
parent, with the luminosity function being much flatter in 
our present model than in that of Cole et al. (2000). At 
Mb — 51og/i « — 13 the difference in amplitude of the B- 
band luminosity functions is about a factor of 4. Similar 
behaviour is seen in the K-band. Compared to the Cole et 
al. model, our new model is in appreciably better agree- 
ment with recent observational determinations in the B and 
K-bands from the 2dFGRS and 2MASS galaxy surveys, by 
Madgwick et al. (2001) and Cole et al. (2001) respectively, 
although in the K-band the predicted slope is still slightly 
too steep. 

Figure ^ also shows the relative importance of the new 
effects we include compared to the Cole et al. (2000) model. 
The dot-dash line shows the effect of including tidal strip- 
ping of satellite galaxies, but no photoionization feedback; 
this is seen to reduce the number of faint galaxies very 
slightly. The dashed line shows the effect of turning on the 
effect of IGM pressure (through the filtering mass) , but not 
the effect of the ionizing background on cooling in halos, 
while the dotted line has the modified cooling in halos turned 
on, but not the filtering mass. Comparing these, we see for 
the photoionization feedback, it is primarily the effects of the 
IGM pressure which suppress the galaxy formation, while 
the reduction in cooling within halos has a smaller effect. 
(Note that these three luminosity functions are calculated 
using the same value of T as for the standard model for the 
purposes of this comparison.) 

In Fig. ^ we show the B-band luminosity function in 
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Figure 10. The bj band galaxy luminosity function at z=0. The 
heavy solid line shows the prediction of our standard model with 
photoionization feedback, and the thin solid line shows the model 
of Cole et al. 2000. The dashed line shows a model with the 
effects of photoionization included but without any feedback from 
supernovae. All model luminosity functions include the effects of 
dust. The symbols show observational data. 



a model where photoionization feedback is included, with 
/esc = 100%, but feedback from supernovae is turned off. 
The star forma tion history for this same model was pre- 
sented in § 1.1.1 . We choose T = 0.95 to match the amplitude 
of the observed luminosity function at L ~ L*. The other 
parameters are the same as in Cole et al. (2000) and in our 
standard model. The low value of T is required because most 
gas which has cooled is locked up into small objects, leav- 
ing little to form bright galaxies. Strictly speaking, a value 
of T < 1 is unphysical, because it requires a negative mass 
in brown dwarfs (defined here as objects with m < O.IM©). 
However, the same results as for T = 0.95 could be obtained 
by small modifications to the IMF at 0.1 < m < IMq, re- 
ducing the mass in low mass stars which anyway contribute 
negligibly to the light from stellar populations. This "no 
SNe feedback" model gives an acceptable match to the ob- 
served luminosity function at the bright end (except possi- 
bly at the highest luminosities) . It predicts a faint-end slope 
which is much steeper than in our standard model, but only 
slightly steeper than the Cole et al. model, which had su- 
pernova feedback but no photoionization feedback, and also 
only slightly steeper than the measurement of Zucca et al. 
(1997). The faint-end slope is still much flatter than in a 
model with no feedback of any type. 

We emphasize that the "no SNe feedback" model we 
have presented here is by no means a "best-flt" model, since 
we have not varied other parameters to achieve a better 
match to the luminosity function, nor have we considered 
other observational constraints as Cole et al. did. (Prelimi- 
nary analysis suggests that a model with only photoioniza- 
tion feedback has difficulties in matching the colours and 
sizes of present-day galaxies). However, the prediction for 
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Figure 9. Galaxy luminosity functions at z=0. The left-hand panel shows the luminosity function in the bj-band while the right-hand 
panel shows that in the K-band. In each panel, the heavy solid line shows the prediction of our standard model with photoionization 
feedback, and the thin solid line shows the model of Cole et al. 2000. Dotted lines show our standard model with the effects of the 
filtering mass switched off, dashed lines show the standard model with photoheating of gas in halos switched off, and dot-dashed lines 
show the standard model with both of these effects switched off. All model luminosity functions include the effects of dust. The symbols 
show observational data. 



the faint-end slope of the luminosity function is expected to 
be fairly robust, so we conclude that if the slope measured 
in the largest and most recent surveys (e.g. Madgwick et al. 
2001) is correct, then photoionization feedback on its own 
does not produce a slope as flat as in the real universe. We 
defer a more detailed study of models without supernova 
feedback to a future paper. 



4.2.2 Tully-Fisher Relation 

Figure |l^ shows the I-band Tully-Fisher relation of galax- 
ies in our model, compared to the observational data of 
Matthewson, Ford & Buchhorn (1992). For constructing 
the model relation, we select galaxies in the same way as in 
Cole et al. (2000), namely we select only star-forming spiral 
galaxies, but also select only those galaxies which have not 
been seriously disrupted by tidal forces (specifically we re- 
move any galaxy which has lost more than 25% of the mass 
of its disk through tidal stripping). These strongly tidally 
disrupted galaxies are unlikely to be recognisable as disks. If 
we do not remove these galaxies, then Tully-Fisher relation 
in our model shows a scatter to very faint magnitudes at low 
circular velocities, but for circular velocities Vc ^lOOkms"^ 
the removal of these galaxies has little effect. 

The figure also shows the model prediction of Cole et al. 
(2000), from which it can be seen that the differences from 
our new standard model are quite small. We also plot lines 
showing the result of switch ing off the effects of photoioniza- 
tion (as described in §4.2.1 and also in the figure caption). It 



can be seen from these curves that tidal stripping makes lit- 
tle difference to the Tully-Fisher relation, while the modified 
cooling in halos has the larger effect at high luminosities and 
the IGM pressure the larger effect at low luminosities. In any 



case, photoionization does not help remove the offset in the 
predicted zero-point of the Tully-Fisher relation relative to 
the observed one. This offset persists to bright magnitudes, 
where photoionization has little effect on galaxy formation. 



4.2.3 Further Properties of the Model at z = 

We now briefiy consider the effects of photoionization feed- 
back on some other predicted properties of galaxies at z = 0. 
We consider the same properties as were compared with ob- 
servational data in Cole et al. (2000). Most of these com- 
parisons concerned fairly luminous galaxies, for which the 
properties in our new model are very similar to those of the 
Cole et al. fiducial model, so we just summarize the main 
results here. 

Cole et al. (2000) computed the distribution of disk 
scale-lengths of spiral galaxies at different luminosities, and 
compared to the observational data of de Jong & Lacey 
(2000), finding good agreement in the magnitude range they 
considered, — 19 > Mi — 5 log h > —22. Our model produces 
almost identical results, as may be expected for these bright 
galaxies. 

Table |l| compares the fractions of S, SO and E galaxies 
in our model brighter than Mb — 51og/i = —19.5 (i.e. L*) 
with Cole et al. (2000) and with observational data. We 
assign morphological types to our model galaxies based on 
their bulge-to-total luminosity ratio in the B-band, B/Tb, 
(including dust extinction). Galaxies having B/Tb < 0.4 are 
classed as S, those with B/Tb > 0.6 are classed as E, and 
those in between are classed as SO. Our model produces a 
slightly higher fraction of spirals than did that of Cole et al. 
(2000), a consequence of the more detailed calculation of 
merger times adopted here. This is in slightly better agree- 
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Figure 11. The TuUy-Fisher relation in the I-band at z = 0. 
The lines show the predicted median relation, with error bars 
indicating the 10% and 90% intervals of the distribution, while 
the filled dots show observational data of Matthcwsoii, Ford & 
Buchhorn (1992). The heavy solid line shows the standard model 
of this paper, while the thin solid line shows the prediction from 
Cole et al. (2000) . Only star-forming spiral galaxies axe included 
in the model relation, with magnitudes corrected to their face- 
on value including the effects of dust. The velocities plotted arc 
the circular velocity at the half-mass radius of the galaxy disc. 
Dotted lines show our standard model with the effects of the 
filtering mass switched off, dashed lines show the standard model 
with photoheating of gas in halos switched off, and dot-daslicd 
lines show the standard model with both of these effects switched 
off. 
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Cole et al. (2000) 


61 


08 


31 


Loveday (1996) 


67 


20 


13 



Table 1 . The morphological mix of galaxies brighter than Mb — 
5 log h = —19.5 from this work and from the model of Cole et al. 
(2000). Also shown is the morphological mix in the APM Bright 
Galaxy Catalogue (which is apparent magnitude limited) from 
Loveday (1996). 

merit with the observational data, but given the crude way 
in which morphological types are assigned in the models, 
these differences should not be over-emphasized. 

The cold gas content of L* spiral and irregular galaxies 
considered by Cole et al. (2000) is unchanged in our model, 
as are the metallicities of gas and stars in these galaxies. 
However, the metallicity of gas in spirals and irregulars does 
show a somewhat steeper dependence on luminosity than in 
Cole et al. (2000), resulting in slightly better agreement 
with observational data. This difference arises because of 
the effect of the filtering mass. In halos only slightly more 
massive than the filtering mass, there can have been no en- 



richment of gas in smaller halos in the merging hierarchy (as 
these halos do not accrete gas). The faint central galaxies of 
these low mass halos are therefore accreting relatively metal 
poor gas compared to those in the Cole et al. (2000) model, 
resulting in lower gas metallicities. The metallicity of stars 
in elliptical galaxies is also changed, but in a different way. 
Bright ellipticals are the same in our model as in that of 
Cole et al. (2000), but faint ones on average have some- 
what higher metallicities than in Cole et al., which worsens 
the agreement with the observational data. Here, the main 
effect is that, with photoionization switched on, a galgixy of a 
given luminosity tends to be found in a halo with higher cir- 
cular velocity (since the filtering mass reduces the amount of 
gas able to accrete into each halo). The higher circular veloc- 
ity implies a deeper potential well, which makes the galaxy 
better at retaining metals (i.e. fewer are lost in the winds 
associated with supernovae feedback), increasing the effec- 
tive yield and raising the metallicities of the low-luminosity 
ellipticals. Pre-processing of gas in lower mass halos is not 
so important for the ellipticals, since the gas is processed 
right up to the effective yield very quickly in the burst of 
star formation which makes the elliptical. 

In conclusion, the largest differences in galaxy proper- 
ties between our new model and that of Cole et al. (2000) 
occur for low luminosity galaxies, the differences beginning 
to be noticeable at around 1 magnitude faintwards of L*. 
The most important difference is a flattening of the faint 
end of the galaxy luminosity function. This slope is in rea- 
sonable agreement with the latest observational estimates 
from the 2dFGRS and 2MASS galaxy surveys. 

4-2.4 Properties of Galaxies at High Redshifts 

Semi-analytic models have been used extensively to investi- 
gate the populations of galaxies seen at high redshifts, such 
as Lymaii-break galaxies (Baugli et al. 1998; Governato et 
al. 1998; Somerville, Primack & Faber 2001). We find that 
photoionization has very little effect on the properties of 
Lyman-break galaxies at 2: = 3, for the range of luminosi- 
ties that is currently observed, because the filtering mass is 
well below the typical mass halo in which these galaxies re- 
side. We defer a more detailed consideration of high-redshift 
galaxies to a future paper. 



5 DISCUSSION 

We have presented a coupled model for evolution of the ion- 
ization state and thermal properties of the IGM and the 
formation of galaxies. The IGM is photoionized by radiation 
from stars in galaxies and from quasars, and the photoioniz- 
ing background in turn exerts a negative feedback effect on 
further galaxy formation. This photoionization feedback op- 
erates in two ways, by heating the IGM, and so by the effects 
of gas pressure reducing the amount of gas which collapses 
into halos, and by ionizing and heating gas within halos, and 
so reducing the amount of gas able to cool to form galaxies. 
The evolution of the ionizing luminosity of the galaxy pop- 
ulation is calculated sclf-consistently with the effects of this 
photoionization feedback. 

We calculate the formation of galaxies within the CDM 
model, by adapting the semi-analytic galaxy formation 
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model of Cole et al. (2000), modified to include the pho- 
toionization feedback effects described above. This is cou- 
pled to a simple model for the evolution of a clumpy IGM, 
which, given the evolution of the ionizing emissivity of galax- 
ies and quasars as an input, predicts the evolution of the 
mean ionized fractions of hydrogen and helium, the volume- 
averaged temperature of the IGM, and the ionizing back- 
ground. We have tested the IGM model against the results 
from numerical simulations of the IGM, and find that the 
predictions for global properties agree reasonably well. In 
particular, we find that our simple IGM model accurately 
predicts the evolution of the characteristic halo mass below 
which accretion of baryonic matter is strongly suppressed, 
which is the most important quantity in our later study of 
galaxy formation. 

In order to more accurately predict the properties of 
satellite galaxies within larger dark matter halos, we have 
also improved the Cole et al. (2000) semi-analytical model 
by incorporating a detailed treatment of the dynamics of 
satellites, including the effects of dynamical friction, tidal 
stripping and heating by tidal shocks. We have compared 
this model in the pure dark matter case with the results from 
high-resolution N-body simulations on the amount of sub- 
structure in dark halos, and find good agreement. This im- 
proved model predicts merging timescales for galaxies that 
on average are comparable to those from the simple esti- 
mates used in previous work, although some satellites have 
their dark halos heavily stripped by tidal forces, and these 
have much longer merging timescales as a consequence of the 
weaker dynamical friction force resulting from the reduced 
satellite mass. 

A significant uncertain parameter in our photoioriiza- 
tion model is /esc, the fraction of ionizing photons from 
stars able to escape from galaxies. In our model, we need 
to assume /esc = 100% in order to produce the emissiv- 
ity of ionizing photons at 2 = 3 inferred observationally by 
Steidel, Pettini & Adelberger (2001). The model then pre- 
dicts reionization of Hi at 2 « 8 and reionization of Hell 
at a w 4. Each reionization event is preceded by an episode 
of reheating. However, this model also produces an ionizing 
background at 2: < 4 which is higher than observational es- 
timates based on the proximity effect. An escape fraction of 
10% gives much better agreement with observational data on 
the ionizing background, but produces an uncomfortably low 
redshift of reionization. In any case, the choice of /esc does 
not change our conclusions about the properties of galaxies 
at 2 = 0. 

Applying our model to the evolution of the galaxy pop- 
ulation, we find the following results: 

(i) The global star formation rate in our model is sup- 
pressed slightly after each episode of reheating due to reion- 
ization of Hi and Hell. The suppression is quite small, with 
reductions of no more than 25% compared to a model with 
no reionization. By z — 0, the star formation rate has re- 
covered to the level predicted by our model with no reion- 
ization, as by then most stax formation is occurring in halos 
well above the masses and temperatures affected by pho- 
toionization feedback. 

(ii) Galaxies brighter than L* are mostly unaffected by 
photoionization. Faintwards of L*, photoionization becomes 
progressively more important, reducing the abundance of 



galaxies of given luminosity. Keeping the same prescrip- 
tion for supernovae feedback as used by Cole et al. (2000), 
we find that including photoionization feedback produces a 
much better fit to recent determinations of the faint end of 
the galaxy luminosity function at 2 = (e.g. Madgwick et al. 
2001; Cole et al. 2001). Most of the effect is due to the in- 
ability of hot IGM gas to accrete into low-mass dark matter 
halos, but heating of gas in halos by the ionizing background 
and tidal limitation of satellite galaxies also play a role. 

(iii) Preliminary analysis of a model with no feedback 
from supernovae, but including the effects of photoioniza- 
tion indicates that such a model can produce a luminosity 
function with faint end slope almost as flat as some observa- 
tional estimates, and significantly fiatter than a model with- 
out supernovae feedback or photoionization. Further work is 
needed to determine if such a model can be made consistent 
with other observational data. 

(iv) Other properties of bright galaxies at the present day 
(e.g. sizes, TuUy-Fisher relation, metallicities) are unaffected 
by photoionization. For faint galaxies, we find differences 
in the Tully-Fisher relation and in metallicities which are 
readily understood. 

(v) Photoionization has little effect on the predicted prop- 
erties of Lyman-break galaxies, over the range of redshifts 
and luminosities for which they are actually observed. These 
galaxies at a = 3 typically live in halos significantly more 
massive than that at which photoionization feedback be- 
comes important, so their properties are insensitive to the 
reionization history. 

If the Universe was reionized through photoionization 
(and no convincing alternative has been proposed), then the 
mechanisms inhibiting galaxy formation which we have ex- 
amined in this paper must operate. As such, no model of 
galaxy formation is complete without their inclusion. Al- 
though we have shown that the properties of bright galaxies 
are almost entirely unaffected, the properties of faint galax- 
ies are strongly infiucnce by photoionization. The methods 
described in this paper provide a flexible and computation- 
ally efficient way to assess the impact of photoionization on 
galaxy formation, and allow us to make definite predictions 
for the properties of faint galaxies. 

As we have shown, photoionization feedback has the 
greatest effect on faint galaxies residing in low mass dark 
matter halos. As such, it will undoubtedly have important 
implications for predictions about the population of satellite 
galaxies found in the Local Group. In the second paper in 
this series, we will explore in detail the properties of these 
galaxies. 
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